cornelius.senf@geo.hu-berlin.de

The 4 R Sessions:

  1. Basic concept of programming, data types, vectors and factors

  2. Matrices and data.frames, data import/export, graphics and visualization

  3. Tidy data, data wrangling

  4. Today: Spatial data in R

Raster data in R

Raster data in R

library(raster)

Raster data in R

ras <- raster(ncol = 100, nrow = 100, xmn = 0, xmx = 100, ymn = 0, ymx = 100)

ras
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Raster data in R

m <- rnorm(n = ncell(ras), mean = 0, sd = 1)

values(ras) <- m

plot(ras)

Raster data in R

Raster data in R

ras2 <- ras + 10

ras3 <- ras * ras2

ras4 <- ras^2

ras5 <- (ras - cellStats(ras, stat = mean)) / 
                                cellStats(ras, stat = sd)

Raster data in R

Raster data in R

ras.na <- ras

ras.na[ras.na > 0.5] <- NA

Raster data in R

Raster data in R

k <- matrix(c(  0, 0.5, 1,
              0.5, 1.5, 2,
              1.5, 6.0, 3), 
            ncol = 3, byrow = TRUE)
k
##      [,1] [,2] [,3]
## [1,]  0.0  0.5    1
## [2,]  0.5  1.5    2
## [3,]  1.5  6.0    3
ras.reclass <- reclassify(x = ras, rcl = k)

Raster data in R

Raster data in R

ras[10, 20]
##           
## 0.7534085
ras[10, 20] <- 10

Raster data in R

Raster data in R

ras.stack <- stack(ras, ras2)

ras.stack
## class       : RasterStack 
## dimensions  : 100, 100, 10000, 2  (nrow, ncol, ncell, nlayers)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## names       :   layer.1,   layer.2 
## min values  : -3.592167,  6.407833 
## max values  :  10.00000,  13.60653

Raster data in R

layer1 <- subset(x = ras.stack, subset = 1)

layer1
## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer.1 
## values      : -3.592167, 10  (min, max)

Raster data in R

data <- as.data.frame(ras.stack)

head(data)
##        layer.1   layer.2
## 1 -0.531517486  9.468483
## 2 -0.183922864  9.816077
## 3  0.003593095 10.003593
## 4 -0.989860242  9.010140
## 5 -3.592166533  6.407833
## 6 -1.795983887  8.204016

Raster data in R

elevation <- raster("Data/elevation.envi")

elevation
## class       : RasterLayer 
## dimensions  : 709, 1306, 925954  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 16.90833, 27.79167, 44.33333, 50.24167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/corneliussenf/Dropbox/Teching/Quantitative Methods/Week-13/Data/elevation.envi 
## names       : elevation

Raster data in R

projection(elevation)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Raster data in R

plot(elevation)

Raster data in R

Raster data in R

# use stack() or brick() for multi-layer raster files
climate <- stack("Data/climate.envi")

# get the number of raster layers
nlayers(climate)
## [1] 4
# get the raster layer names
names(climate)
## [1] "X.Tmax."  "X.Tmin."  "X.Tmean." "X.Prec."
# rename the raster layers
names(climate) <- c("Tmax", "Tmin", "Tmean", "Prec")

Export raster data

writeRaster(elevation, 
            filename= "Data/elev.envi",
            overwrite = TRUE)

Vector data

Vector data

library(sp)

Vector data: simple random sample

coords_srs <- data.frame(x = runif(100, 17, 27), 
                         y = runif(100, 45, 50))

pts_srs <- SpatialPoints(coords = coords_srs)
pts_srs
## class       : SpatialPoints 
## features    : 100 
## extent      : 17.04253, 26.99759, 45.01899, 49.98391  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

Vector data: Assign Coordinate system

proj <- projection(elevation)

CRS(proj)
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
pts_srs <- SpatialPoints(coords = coords_srs,
                    proj4string = CRS(proj))

Vector data

plot(elevation)
plot(pts_srs, add = TRUE)

Advanced plotting

library(mapview)

mapview(climate)

Advanced plotting

Vector data

df.climate <- extract(climate, pts_srs, df = TRUE)

class(df.climate)
## [1] "data.frame"
head(df.climate)
##   ID     Tmax     Tmin     Tmean     Prec
## 1  1 161.3333 60.66667 110.66666 47.83333
## 2  2 137.1667 38.58333  87.58334 47.33333
## 3  3 118.1667 38.00000  77.83334 56.33333
## 4  4 154.0000 53.91667 103.58334 46.75000
## 5  5 116.8333 28.33333  72.41666 65.33334
## 6  6 147.9167 40.50000  93.91666 48.50000

Vector data

spdf.climate <- extract(climate, pts_srs, sp = TRUE)

spdf.climate
## class       : SpatialPointsDataFrame 
## features    : 100 
## extent      : 17.04253, 26.99759, 45.01899, 49.98391  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 4
## names       :            Tmax,             Tmin,            Tmean,             Prec 
## min values  :           83.75,             -2.5, 40.3333320617676, 43.4166679382324 
## max values  : 164.16667175293, 65.6666641235352,           114.25, 86.4166641235352

Vector data

You can acess SpatialPointsDataFrame the same way as data.frame!

head(spdf.climate)
##       Tmax     Tmin     Tmean     Prec
## 1 161.3333 60.66667 110.66666 47.83333
## 2 137.1667 38.58333  87.58334 47.33333
## 3 118.1667 38.00000  77.83334 56.33333
## 4 154.0000 53.91667 103.58334 46.75000
## 5 116.8333 28.33333  72.41666 65.33334
## 6 147.9167 40.50000  93.91666 48.50000
spdf.climate$Tmax[1]
## [1] 161.3333

Reading vector data

library(rgdal)

Reading vector data

pts.sample <- readOGR('Data', 'pts_sample')
## OGR data source with driver: ESRI Shapefile 
## Source: "Data", layer: "pts_sample"
## with 142 features
## It has 1 fields
## Integer64 fields read as strings:  id
pts.sample
## class       : SpatialPointsDataFrame 
## features    : 142 
## extent      : 17.80209, 27.04407, 44.64419, 49.84834  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 1
## names       : id 
## min values  :  0 
## max values  : 99

Writing vector data

writeOGR(spdf.climate, "Data", "pts_climate", driver = "ESRI Shapefile", overwrite_layer = TRUE)